Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 30 October 2009 (MN KTeX style file v2.2) 



Chaotic Transport and Chronology of Complex Asteroid Families 



o 
o 

> 

o 
O 

o 

m 

o 



> 

in 

00 

d 

On 
O 



B. Novakovic^*, K. Tsiganis^* and Z. Knezevic^* 

^ Astronomical Observatory, Volgina 7, 110 60 Belgrade 38, Serbia 

^ Section of Astrophysics, Astronomy & Mechanics, Department of Physics, Aristotle University of Thessaloniki, GR 54 124 Thessaloniki, Greece 



30 October 2009 



ABSTRACT 

We present a transport model that describes the orbital diffusion of asteroids in chaotic regions 
of the 3-D space of proper elements. Our goal is to use a simple random-walk model to study 
the evolution and derive accurate age estimates for dynamically complex asteroid families. To 
this purpose, we first compute local diffusion coefficients, which characterize chaotic diffu- 
sion in proper eccentricity (ep) and inclination {Ip), in a selected phase-space region. Then, a 
Monte-Carlo-type code is constructed and used to track the evolution of random walkers (i.e. 
asteroids), by coupling diffusion in (epjp) with a drift in proper semi-major axis (a^) induced 
by the YarkovskyA'ORP thermal effects. We validate our model by applying it to the family 
of (490) Veritas, for which we recover previous estimates of its age (^ 8.3 Myr). Moreover, 
we show that the spreading of chaotic family members in proper elements space is well repro- 
duced in our randomk-walk simulations. Finally, we apply our model to the family of (3556) 
Lixiaohua, which is much older than Veritas and thus much more affected by thermal forces. 
We find the age of the Lixiaohua family to be 155 ± 36 Myr. 

Key words: celestial mechanics, minor planets, asteroids, methods: numerical 



1 INTRODUCTION 



As first noted by iHiravamal ( Il918h . asteroids can form prominent 
' groupings in the space of orbital elements. These groups, nowa- 
days well-known as asteroid families, are believed to have resulted 
from catastrophic collisions among asteroids, which lead to the 
ejection of fragments into nearby heliocentric orbits, with relative 
■ velocities much lower than their orbital speeds. To date, several 
' tens of familie s have been discovered across the whole asteroi d 
' Main Beh (e.g. lSendiova & ZapDalil2002l : iNesvornv et al.ll2006h. 
Also, families have been i dentified among the Trojans jMilanil 
Il993l : iBeauge & Roidl200lh . and most recently, pro posed to ex- 
ist in the Transneptunian region l lBrown etal J 120071) . Studies of 
asteroid families are very important for planetary science. Fami- 
lies can be used, e.g. to understand the coUisional history of the 
asteroid Main Belt (Bottke et al. 2005,) , the outcomes of disrup- 
tion events over a size range inaccessible to laboratory experiments 
(e.g. lMichel et alj2003l : lDurda et^l2007h . to understand the min- 
eralogical structure of their parent bodies (e.g. lCellino et alj|2002t) 
and t he effects of related dust "showers" on the Earth jparlev et al.l 
I2OOQ) . Obtaining the relevant information is, however, not easy. 
One of the main complications arises from the fact that the age of 
a family is, in general, unknown. Thus, accurate dating of asteroid 
families is an important issue in the asteroid science. 

A number of age-determination methods have been proposed 
so far. Probably the most accurate procedure, particularly suited for 
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young families (i.e. age < 10 Myr), is to integrate the orbits of the 
family members backwards in time, until the orbital orientation an- 
gles cluster around some value. As such a conjunction of the orbital 
elements can occur only immediately after the disruption of the par- 
ent body, the time of conjunction indi cates the format i on time. The 
method was successfully applied by INesvomv et alj ( I2OO2I . l2003h 
to estimate the ages of the Karin cluster (5.8±0.2 Myr) and of the 
Veritas family (8.3±0.5 Myr). This method is however limited to 
groups of objects residing on regular orbits. 

For older families (i.e. age > 100 Myr), one can make use 
of the fact that asteroids slowly spread in semi-major axis due to 
the ac tion of Yarkovsky thermal forces (Farinella & Vokrouhlick^ 
1 19991) . As small bodies drift faster than large bodies, the distribu- 
tion of family members in the {ap,H) plane - where ap is the 
proper semi-major axis and H is the absol ute magnitude - can be 
used as a clock. That method was used bv lNesvomv et al.l j2005l) 
to estimate ages of many asteroid families. In these estimations the 
initial sizes of the families were neglected, so that this methodol- 
ogy can overestimate the real age by a factor of as much as 1.5 -2. 
An improved version of this method, which accounts for the initial 
ejection velocity field and the action of YORP thermal torq ues, has 
been su cces sfully applied to several fami lies by Vokrouh l ickv et al.l 
ilOOei) and lVokrouhlickvetal.1 j2006bh . A gain, it is not straight- 
forward to apply this method to families located in the chaotic re- 
gions of the asteroid belt. 

iMilani & Farinellal jl994) suggested that asteroid families, 
which reside in chaotic zones, can be approximately dated by 
chaotic chronology. This method is based on the fact that the age 
of the family cannot be greater than the time needed for its most 
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chaotic members to escape from the family region. In its original 
form, this method provides only an upper bound for the age. Re- 
centlv. lrsiganis et aD ( l2007h introduced an improved version of this 
metho d, based on a s t atistic al description of transport in the phase 
space. iTsiganis et al.l jlOOTh successfully applied it to the family 
of (490) Veritas, finding an age of 8.7±1.7 Myr, v yhich is statisti- 
cally t he same as that of 8.3±0.5 Myr, obtained bv lNesvomv et al.l 
(l2003fl Despite these improvements, the chaotic chronology still 
suffered from two important limitations. It did not account for the 
variations in diffusion in different parts of a chaotic zone, which can 
significantly alter the distribution of family members (i.e. the shape 
of the family). Moreover, it did not account for Yarkovsky/YORP 
effects, thus being inadequate for the study of older families. 

In this paper we extend the chaotic chronology method, by 
constructing a more advanced transport model, which alleviates 
the above limitations. We first use the Veritas family as bench- 
mark, since its age can be considered well-defined. Local diffu- 
sion coefficients are numerically computed, throughout the region 
of proper elements occupied by the family. These local coeffi- 
cients characterize the efficiency of chaotic transport at different 
locations within the considered zone. A Monte-Carlo-type model 
is then constructed, in analogy to the one used by Tsiganis et al. 
( |2007[) . The novelty of the present model is that it assumes vari- 
able transport coefficients, as well as a drift in semi-major axis due 
to Yarkovsky/YORP effects, although the latter is ignored when 
studying the Veritas family. Applying our model to Veritas, we find 
that both (a) the shape of its chaotic component and (b) its age 
are correctly recovered. We then apply our model to the family of 
(3556) Lixiaohua, another outer-belt family but much older than 
Veritas and hence much more affected by the Yarkovsky/YORP 
thermal effects. We find the age of the Lixiaohua family to be 
~ 155 Myr. 

We note that, depending on the variability of diffusion coeffi- 
cients in the considered region of proper elements, this new trans- 
port model ca n be computationally much more expensive than the 
one applied in^ Tsiganis et al.l ( |2007|) . This is because, if the values 
of the diffusion coefficients vary a lot across the considered region, 
one would have to calculate them in many different points. How- 
ever, even so, this computation needs to be performed only once. 
Then, the random-walk model can be used to perform multiple runs 
at very low cost, e.g. to test different hypotheses about the original 
ejection velocities field or about the physical properties of the as- 
teroids. On the other hand, for "smooth" diffusion regions in which 
the coefficients only change by a factor of 2-3 across the consid- 
ered domain, the model can be simplified. In such regions, the age 
of a family can be accurately determined even by assuming an av- 
erage (i.e. constant over the entire region) diffusion coefficient, as 
we show in Section 3. 



2 THE MODEL 

Our study begins by selecting the target phase-space region. This is 
done by identifying the members of an asteroid family crossed by 
resonances, from a catalog of numbered asteroids. Apart from the 
largest Hirayama families, for the other smaller and more compact 
ones, in the current catalog one typically finds up to several hun- 
dred members. Thus the chaotic component of the family consists 

^ Of course, Tsiganis et al. (2007) and Nevorny et al. (2003) used a chaotic 
and a regular subsets of the family, respectively. 



of a few tens to a few hundreds of asteroids. Although this may be 
adequate to compute the average values of the diffusion coefficients 
(as in Tsiganis et al. 2007), a detailed investigation of the local dif- 
fusion characteristics requires a much larger sample of bodies. The 
latter can be obtained by adding in the fictitious bodies, selected in 
such a way that they occupy the same region of the proper elements 
space as the real family members. Since we wish to study just these 
local diffusion properties and the effect of the use of variable coef- 
ficients in our chronology method, we are going to follow here this 
strategy. 

The 3-D space of proper elements, occupied by the selected 
family, is divided into a number of cells. Then, in each cell, the dif- 
fusion coefficients are calculated for both relevant action variables, 
namely Ji ~ 5 \J TTj^i ~ \\J TTi ™^ ^'^J denotes 

Jupiter's semi-major axis, the proper eccentricity and Ip the 
proper inclination of the asteroid). This is done by calculating the 
time evolution of the mean squared displacement {(AJi)^) (i=l,2) 
in each action, the average taken over the set of bodies (real or fic- 
titious) that reside in this cell. The diffusion coefficient is then de- 
fined as the least-squares-fit slope of the ((AJi) ^) (t) c urve, while 
the formal error is computed as in lTsiganis et al] faoTl) . 

The simulation of the spreading of family members in the 
space of proper actions and the determination of the age of the fam- 
ily is done using a Markov Chain Monte Carlo (MCMC) technique 
(e.g. Gentle 2003; Berg 2004). At each step in the simulation the 
random walkers can change their position in all three directions, i.e. 
the proper semi-major axis ap and the two actions Ji and J2. Al- 
though no macroscopic diffusion occurs in proper semi-major axis, 
the random walker can change its ap value due to the Yarkovsky 
effect, while the changes in Ji and J2 are controlled by the local 
values of the diffusion coefficients. In the case of normal diffusion 
the transport properties in action space are determined by the solu- 
tion o f the Fokker-Planck equation (see lLichtenberg & LiebermanI 
(l983)). The MCMC method is in fact equivalent to solving a 
discretized 2-D Fokker-Planck equation with variable coefficients, 
combined here with a 1-D equation for the Yarkovsky-induced dis- 
placement in ap. The latter acts as a drift term, contributing to the 
variability of diffusion in Ji and J2 . 

The rate of change of ap due to the Yarkovsky thermal 
force, is given by the following equation (e.g. IVokrouhlickvlll999l: 
iFarinella & VokrouhUck^ll999r) : 

— = fci cos 7 + ^2 sin^ 7 (1) 
dt 

where the coefficients k\ and ^2 depend on parameters that de- 
scribe physical and thermal characteristics of the asteroid and 7 de- 
notes the obliquity of the body's spin axis. For km-sized asteroids, 
the drift rate is inversely proportional to their radius. This simpli- 
fied Yarkovsky model assumes that the asteroid follows a circular 
orbit, and thus linear analysis can be used to describe heat diffusion 
across the asteroid's surface. The obliquity of the spin axis and the 
angular velocity of rotation (uj) of the asteroid are subject to ther- 
mal torques (YORP) that change their values with time, according 
to the following equations: 

(e.g. IVokrouhUckv & CaDekll2002l : ICapek & Vokrouhlickvl l2004l) . 
where the functions / and g describe the mean strength of the 
YORP to rque and depend on the aster oid's surface thermal con- 
ductivity dCapek & VokrouhUckvll2004h . 

The length of the jump in ap that a random-walker under- 
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takes at each time-step dt in the MCMC simulation, is determined 
by equations (l)-(2), in their discretized form. Of course, a set 
of values of the physical parameters must be assigned to each 
body. As the majority of the Veritas family member s are of C- 
type jPi Martino et alJI 19971 : iMothe-Diniz et al.ll2005h. while the 
Lixiaohua family members seem to be C/X-tvpe jLazzaro et al.l 
|2004 I Nesvomv et al.l |2005[), the follo wing values for these pa- 
rameters i Broz et al. 20051 : lBro3l2006b are adopted: thermal con- 
ductivity K = 0.01 — 0.5 [W(m K)~^], specific heat capacity 
C = 1000 [J(Kkg)^^], and t he same value f or sur face and bulk 
density g — 1500 [kgm' "^]. In lTedesco et"aLl ( l2002h . the geomet- 
ric albedos of several Veritas and Lixiaohua family members 
are listed, yielding a mean p„ = 0.068 ± 0.018 for Veritas and 
Pv = 0.049 ± 0.027 for Lixiaohua. The rotation period, P, is cho- 
sen randomly from a Gaussian distribution peaked at P = 8 h, 
while the distribution of initial obliquities, 7, is assumed to be uni- 
fornfl To assign the appropriate values of absolute magnitude H 
to each body, we need to have an estimate of the cumulative distri- 
bution A^(< H) of family members. A power-law approximation 
is used (e.g. lVokrouhlickv et al.l2006bl) 



iV(< H) oc 10 



(3) 



where /3 depends on the considered interval for H; e.g. for the Ver- 
itas family, we find (5 = 0.74 ± 0.03 for H G[l 1.5,13.5] and 
13 = 0.23 ± 0.03 for H £[13.5,15.5]. Having the values of H 
and Pn, the radius _R of a body can be estimated, using the relation 





(e.g. lCarruba et alj|2003h 



R (km) = 1329 



10— 



(4) 



At each time-step in the MCMC simulation, a random- walker 
suffers a jump in Ji and J2, whose length is given by A J; 
=yL\/ D{Ji)dt/2 (i = 1, 2), where u, is a random number from a 
Gaussian distribution jTsiganis et al.l2007h . Since the values of the 
diffusion coefficients D{Ji) vary in space, the maximum allowable 
jump, for a given dt, changes from cell to cell. In our simulations, 
the values of D{Ji) and -D( J2) used for each body, are given by: 



dn 



dr 



-Df. 



(5) 



dh+dn ' dL+ dR 

where and da denote the distances of the random- walker from 
the two nearest nodes (left and right) and Dl and Dr denote the 
corresponding values of the diffusion coefficients at these nodes|f] 
For a correct determination of the age of the family, the ran- 
dom walkers have to be placed initially in a region, whose size is 
as close as possible to the size that the real family members oc- 
cupied, immediately after the family-forming event. This is in fact 
a source of uncertainty for our model. In our calculations we as- 
sumed the initial spread of the family in (ap,ep) and (ap,sinJp) 
to be accurately represented by a Gaussian equivelocity ellipse (see 
Morbidelli et al. 1995), computed such that (i) the spread in ap 
of the whole family and (ii) the spread in Cp and sin Ip of family 
members that follow regular orbits is well reproduced. 



^ According to lPaolicchi et alj il99d) a size-rotation relation suggests that 
sma ller fragments a re rota ting somewhat faster than the larger ones (see 
also lDonnisonl2003h . Also. lKrvszczvnska et alj i2007h claim that the poles 
are not isotropically distributed, as general theoretical considerations may 
predict. These facts are not considered here, but could become important 
when studying very old families. 

^ A geometric mean D{Ji)=^Di^DR could be used instead of an arith- 
metic one; we actually found negligible differences. 



3 THE RESULTS 

In this section we use our model to study the evolution of the 
chaotic component of two outer-belt asteroid families: (490) Ver- 
itas and (3556) Lixiaohua. In both cases, a number of mean motion 
resonances (MMR) cut-through the family, such that a significant 
fraction of members follow chaotic trajectories. On the other hand, 
their ages differ significantly, according to previous estimates. In 
this respect, the Yarkovsky effect can be neglected in the study of 
Veritas, but not in the study of Lixiaohua. 

We begin by performing an extensive study of the local dif- 
fusion properties in the chaotic region of the Veritas family. Then, 
the MCMC model is used to simulate the evolution of the chaotic 
members and to derive an estimate of the age of the family. The re- 
sults are compared to the ones given by the model of Tsiga nis et al.l 
( 2007). Finally, we apply the MCMC model to the Lixiaohua fam- 
ily and derive estimates of its age, for different values of the 
Yarkovsky-related physical parameters. 



3.1 The Veritas family 

The Veritas family is a comparatively small and compact outer- 
belt family, spectroscopically different from the background popu- 
lation of asteroids. In terms of dynamics, it occupies a very inter- 
esting and complex region, crossed by several mean motion reso- 
nances. Ap plication of the Hierarchical Clustering Method (HCM) 
IZappala et al. (1995) to the AstDys catalog of synthetic proper ele- 
ments (numbered asteroids http://hamilton.unipi.it/astdys as of De- 
cember 2007), yield s 409 family members , for a velocity cut-off of 
Vc = 40 m s~^ as in lTsiganis et al.l j2007h . 

Although the family appears now to extend beyond ap — 
3.18 AU (see Fig.[I), the main dynamical groups remain practically 
the same (see Tsiganis et al. 2007, for a detailed description of the 
groups). Since the scope of this paper is to present a refined trans- 
port model, we will briefly describe here only the main relevant 
features, referring to a forthcoming paper for a renewed analysis of 
the Veritas family itself. 

The main chaotic zone, where appreciable diffusion in proper 
elements is observed, is located around 3.174 AU (Fig. 1) and is 
associated with the action of the (5,-2,-2) three-body mean motion 
resonance (MMR); see Fig. [2] for the typical short-term evolution 
of such a resonant asteroid. The family members that reside in this 
resonance can disperse over the observed range in Cp and sin Ip on 
a ~ 10 Myr time-scale. This is exactly the group of bodies (group 
A) that was used by Tsiganis et al. (2007), to compute the age of 
Veritas. 



3.1.1 The local diffusion coefficients 

As the number of bodies in group A is small, we need to generate a 
uniform distribution of fictitious bodies, in order to compute local 
diffusion coefficients across the observed range in (ep, sin Ip). For 
this reason we start by selecting ~ 25, 000 initial conditions (fic- 
titious bodies), covering the same region as the real Veritas family 
members, in the space of osculating elements. We note that the ac- 
tual number of bodies used in the calculations of the coefficients 
is much smaller than that (see below). The orbits of the fictitious 
bodies are integrated for a time-span of 10 Myr, using the 0RBIT9 
integrator (version 9e), in a model that includes the four major plan- 
ets (Jupiter to Neptune) as perturbing bodies. The indirect effect of 
the inner planets is accounted for by applying a barycentric correc- 
tion to the initial conditions. This model is adequate for studying 
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Figure 1. Distribution of tlie real Veritas family members in the (ap,ep) 
and (ap,sinlp) planes. The superimposed ellipses repre sent equivelocity 
curve s, computed according to the equations of Gauss (e.g. lMorbidelli et aU 
for a velocity of 40m s ^, true anomaly / = 30° and argument 
of pericentre w = 150°. The vertical dashed lines represent approximate 
borders of the main three-body MMRs, as indicated by the corresponding 
labels. 

outer-belt asteroids. Note that the integration time used here is in 
fact longer than the known age of the Veritas family. This is done in 
order to study the convergence of the computation of the diffusion 
coefficients, with respect to the integration time-span. 

For each body mean elements are computed on-line, by ap- 
plying digital filtering, and proper elements ar e subs equently com- 
puted accor ding to the analytical theo ry of iMilan i & Knezevic 
( Il990lll994h . Synthetic proper elements j Knezevic & Milaniii2000l) 
are also calculated, for comparison and control. Since the mapping 
from osculating to proper elements is not linear, the distribution of 
the fictitious bodies in the space of proper elements is not uniform, 
which can complicate the statistics. A smaller sample of ~ 10, 000 
bodies, with practically uniform distribution in proper elements, is 
therefore chosen. Thus, our statistical sample, on which all com- 
putations are based, is in fact ~ 25 times larger than the actual 
population of the family. 

As explained in the previous section, we computed local dif- 
fusion coefficients, by dividing the space occupied by the Veritas 
family in a number of cells. Our preliminary experiments suggested 
that, while a large number of cells is needed to accurately represent 
the dependence of the coefficients on ap, the same is not true for 
Ep and sin Ip, except for the wide chaotic zone of the (5, —2, —2) 
MMR. Thus, we decided to follow the strategy of using a large 
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Figure 2. The time evolution of the mean semi-major axis a (top) and the 
resonant angle (bottom), for a fictitious body inside the (5, -2, -2) three-body 
MMR. Note the correlation between the two quantities. Temporary capture 
at the resonance border occurs when a is at maximum (resp. minimum) and 
a"(5 _2,_2) circulates in a positive (resp. negative) sense. 



number of cells in Up and a small number of cells in and sin Ip, 
except in the (5,-2,-2) region. The efficiency of the computation is 
improved if we use a moving-average technique (i.e. overlapping 
cells), instead of a large number of static cells, because in the latter 
case we would need a significantly larger number of fictitious bod- 
ies. We selected the size of a cell in each dimension as well as an ap- 
propriate step-size, by which we shift the cell through the family, as 
follows: for ap, the cell-size was Aap — 5x 10"** AU and the step- 
size 10"'* AU; for Ji the cell-size was A Ji — 10"'* and the step- 
size 4 X 10~^; finally, for J2, the cell-size was AJ2 = 2 x 10~* 
and the step-size 10~*. Thus, the total number of (overlapping) 
cells used in our computations was 2, 196. 

The time evolution of the mean squared displacement in Ji 
and J2 is shown in Fig. [3] for a representative cell in the (5,-2,-2) 
resonance. The evolution is basically linear in time, as it should be 
for normal diffusion. The slope of the fitted line defines the value 
of the local diffusion coefficient. When performing such computa- 
tions, one needs to know (i) what is the shortest possible integration 
time-span, and (ii) what is the smallest possible number of fictitious 
bodies per cell, for which reliable values of the coefficients can be 
obtained. For several different groups of fictitious bodies (i.e. dif- 
ferent cells), we calculated the diffusion coefficients using differ- 
ent values of the integration time-span, between 1 and 10 Myr. Our 
results suggest, that an integration time of ~ 4 Myr is sufficient 
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Figure 3. The mean squared displacement in Ji and J2 as a function of 
time, for a cell inside the (5,-2,-2) MMR. The evolution is basically linear 
in time, as can be seen from the respective fit given on each graph, with a 
superimposed small amplitude oscillation of a period of ~ 5 Myr. 



to obtain reliable values, as shown in Fig. [4] This saturation time 
is about half the known age of the Veritas family. Hence, for this 
case, computing diffusion coefficients is practically as expensive 
as studying the evolution of the family by long-term integrations. 
However, the saturation time is related to the resonance in ques- 
tion and not to the age of the family, which could be much longer. 
Thus, as a matter of principle, the computational gain can become 
important when dealing with much older families. Resonances of 
similar ord er are characterized by sim ilar Lyapunov and diffusion 
times (see jMurrav & Holmanlll997l) ). and thus similar computa- 
tion time-spans (i.e. a few Myr) should be used, for various reso- 
nances throughout the belt. 

The dependence of the diffusion coefficients on the number 
of bodies considered in each cell was also tested. For a number 
of different cells, we calculated the coefficients, using from 10 to 
100 bodies for the computation of the corresponding averages. Our 
results suggest that at least 50 bodies per cell are needed, for an 
accurate computation. 

The values of the diffusion coefficients in Ji and J2, along 
with their formal errors, are given as functions of in Fig.|5] The 
largest diffusion rate is measured in the (5,-2,-2) MMR, which cuts 
through the family at Op ~ 3.174 AU. Both coefficients increase as 
the center of the resonance is approached, but show local minima 
at flp ~3.174 AU, which is approximately the location of the cen- 
ter. The maximum values are I>(Ji) = (1.28±0.02)x 10"" yr"\ 
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Figure 4. The values of the diffusion coefficients -D(Ji) and D{J2) as a 
function of the integration time span. Each curve on these graphs corre- 
sponds to a different cell. 



local minima are L>(Ji) = (0.73±0.01)x 10"" yr-\ and D{J2) 
= (1.16±0.02)x 10"^'* yr"^- Thi s form of dependence in i s 
in agreement with the results of iNesvomv & Morbidellil ( Il999h . 
where it was shown that the dynamics in this resonance are simi- 
lar to those of a modulated pendulum (see also Morbidelli ( 2002)), 
with an island of regular motion persisting at the center of the res- 
onance. However, the size of the island decreases with decreasing 
Bp, so that the diffusion coefficients may decrease as the resonance 
center is approached, but do not go to zero. 

As also seen in Fig.|5] the (5,-2,-2) is by far the most important 
resonance in the Veritas region, associated with the widest chaotic 
zone. Bodies inside this resona nce exhibi t a com plex behaviour, 
as already noted by Knezevic & Pavlovi^ ( l2002h . Most resonant 
bodies show oscillations in proper semi-major axis around ap = 
3.174 AU, but some are temporarily trapped near the resonance's 
borders (see also Fig.|2]l, at Op = 3. 172 AU or ap = 3. 1755 AU. This 
"stickiness" can be important, as it can affect the diffusion rate. In 
fact, we find that the Op values of these bodies shift towards the 
resonance borders, where slower diffusion rates are also measured. 

The diffusion properties are quite different at the (3,3,-2) (lo- 
cated at ap ^ 3.168 AU) and (7,-7,-2) (at Op « 3.18 AU) MMRs, 
which are of higher order in eccentricity with respect to the (5,-2,- 
2) MMR (see Nesvorny & Morbidelli 1999). The values of D{Ji) 
for the (3,3,-2) resonance (see Fig. |5jt) are also increasing as the 
center of the resonance is approached, but no local minimum is 
seen near the center of the resonance, at least in this resolution. 



and D{J2) = (1.36±0.02)x 10~" yr"\ while the values at the The maximum values are only D{Ji) = (8.30±0.09)x lO"^'' yr" 
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Figure 5. The local diffusion coefficients D(J\) and D{J2) (with their cor- 
responding error-bars) in the Veritas family region, shown here as functions 
of the proper semi-major axis a^. The embedded rectangles are magnifica- 
tions of the graph, in the vicinity of the (3, 3. -2) and (7, -7, 2) MMRs. Note 
that 15 (J2) is practically zero everywhere, except in the (5,-2,-2) region. 



and D( J2) = (0.22±0.03)x 10"^" yr~\ clearly much smaller than 
in the (5,-2,-2) MMR. Note that D{J2) has practically zero value. 
The region of the (7,-7,-2) resonance is even less exciting. D{Ji) 
is almost constant across the resonance, with a very small value 
(4.1±0.06 X 10"^^ yr"^), and D( J2) is practically zero. 

The above results suggest that the (5,-2,-2) MMR is essentially 
the only resonance in the Veritas region characterized by apprecia- 
ble macroscopic diffusion. We now focus on the variation of the dif- 
fusion coefficients with respect to (e^, sin Ip) along this resonance. 
As shown in Fig.|6l the dependence of the diffusion rate on the ini- 
tial values of the actionQ is very complex. The values of D{Ji) 
vary from (0.60±0.01)x 10"" yr'^ to (1.66±0.02)x 10~" yr"^ 



, they vary from (0.63±0.02)x 10" 

-14 „_-l 



yr 



to 



while, for D{J2 

(2.31±0.03)x 10~^^ yi~^ ■ Consequently, chaotic diffusion along 
this resonance can in principle produce asymmetric "tails" in the 
distribution of group- A members. Note, however, that the coeffi- 
cients only vary by a factor of 2 — 3 and that thei r average values 
are essentially the same as in lTsiganis et al.l j2007l) . 

3. 1.2 The MCMC Simulations - Chronology of Veritas 

We can now use the MCMC method to simulate the evolution and 
determine the age of the Veritas family, assuming that all the dy- 
namically distinct groups originated from a single brake-up event. 

* The values of the diffusion coefficients are calculated for a regular grid 
in Ji and J2 and then translated into proper elements space. 
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Figure 6. The local diffusion coefficients D{J\) (top) and £>( J2) (bottom) 
inside the (5,-2,-2) region, shown here in grey-scale as functions of Cp and 
sin /p. The values are averaged over the resonant range in semi-major axis. 



A set of six values (ap,Ji,J2,P,^,H) is assigned to each random 
walker in the simulation. All bodies are initially distributed uni- 
formly inside a region of predefined size in SJi{0) and semi-major 
axes in the range [3.172, 3.176]. The age of the family, (r), is de- 
fined as the time needed for 0.3% of the random walkers to leave 
an ellipse in the (Ji, J2) plane, corresponding to a 3-a confidence 
interval of a 2-D Gaussian distribution. We note that, for Veritas, 
the mobility in semi-major axis due to Yarkovsky is very small and 
practically insignificant for what concerns the estimation of its age, 
since the family is young and distant from the Sun. For older fami- 
lies one should also define appropriate borders in Up. 

Using the values of the coefficients obtained above, and the 
values of a( Ji)=(2.31±0.22)x 10"^, a( J2)=(3.97±0.30)x 10"*, 
calculated from the distribution of the real group A members, we 
simulate the spreading of group A and estimate its age. Of course, 
the model depends on some free parameters: the initial spread of 
the group in (<5Ji(0), (5J2(0)), the time-step, dt, and the number 
of random-walkers, n. Therefore, the dependence of the age, r, on 
these parameters was checked. Uncertainties in the values of the 
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current borders of the group (i.e. the confidence ellipse) and the 
values of the diffusion coefficients were taken into account, when 
calculating the formal error in r. 

Different sets of simulations were performed, the results of 
which are given in Fig.|7ja)-(d). Each "simulation" (i.e. each point 
in a plot) actually consists of 100 different realizations (runs) of 
the MCMC code. In each run, the values of D{Ji) and a(Ji) were 
varied, according to the previously computed distributions of their 
values. The values of the free parameters were the same for all runs 
in a given simulation. 

The first set of simulations was performed in order to check 
how the results depend on the time step, dt. Five simulations were 
made, with dt ranging from 1000 to 5000 yr (Fig.EJi). The stan- 
dard deviation of r is relatively small, suggesting that r is roughly 
independent of dt. According to this set of simulations, the age of 
the family is r = (r) ± At — 8.5 ± 1.3 Myr, where (t) is the 
mean value and At the standard error of the megin. This esti mate 
is in excellent agreem ent with those of iNesvomv et al.l ( |2003|) and 
lTsiganisetal]j2007l) . 

The second set of simulations was performed in order to check 
how the results depend on the number of random walkers, n. As 
shown in Fig.[7]3), r is weakly dependent of n, the variation of the 
mean value of r is slightly larger than in the previous case. The age 
of the family, according to this set, is r = 8.6±1.3 Myr. 

The third group of simulations was performed in order to 
check how the results depend on the assumed initial spread of group 
A in (5J2(0), a parameter which is poorly constrained from the re- 
spective equivelocity ellipse in Fig.[T] We fixed the value of 5Ji (0) 
to 2.3x10"*, which can be considered an upper limit, according 
to Fig. (ijfl- Six simulations were performed, with 5J2(0) ranging 
from 3.5 xlO"'' to 11.0x10"*, and the results are shown in Fis. 
[TJ;. The values of r tend to decrease as SJ2 (0) increases. This is to 
be expected, since increasing the initial spread of the family, while 
targeting for the same final spread, should take a shorter time for a 
given diffusion rate. The results yield r = 8.8±1.1 Myr. 

As a final check, we performed a set of simulations with 
5Ji(0) = 2.3x10"* and (5J2(0) = 11.0x10"*. These values cor- 
respond to equivelocity ellipses that contain almost all regular and 
(3,3,-2)-resonant family members, except for the very low inclina- 
tion bodies (sin Jp < 0.156). Five sets of runs, for five different 
values of dt, were performed (see Fig.|7jl). As in our first set of 
simulations, r is practically independent of dt. On the other hand, 
T turns out to be smaller than in the previous simulations, since the 
assumed values for SJi (0) and 5 J2 (0) are quite large. Even so, we 
find r = 7.6±1.1 Myr, which is still an acceptable value. 

Combining the results of the first three sets of simulations and 
taking into account all uncertainties, we find an age estimate of r = 
8.7±1.2 Myr for the Veritas family . This result is very close to the 
one found bv lTsiganis et al.l ilOoH) . the error though being smaller 
by ~ 30%. 

In addition to the determination of the family's age, we would 
like to know how well the MCMC model reproduces the evolu- 
tion of the spread of group-A bodies, in the (Ji, J2) space. For this 
purpose we compared the evolution of group A for 10 Myr in the 
future, as given (i) by direct numerical integration of the orbits, and 
(ii) by an MCMC simulation with variable diffusion coefficients. 
Figures[Ha)-(b) show the outcome of this comparison. As shown in 



° Assuming that an equivelocity ellipse in (ap,ep), large enough to encom- 
pass both the regular part of the family and the (3,3,-2) bodies, is a better 
constraint. 
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Figure 7. Dependence of r and (t(t) on three free parameters: (a) time step 
dt; (b) number of random walkers n; (c) initial size of the family <5J2(0). 
The bottom panel, 7(d), shows the dependence of r on dt for a run in which 
the maximum values of 5Ji(0) and <5J2(0) were assumed. The values of 
the free parameters that ai'e constant in each group of simulations ai'e indi- 
cated on top of each plot. The horizontal lines and dashed areas denote the 
age of the Veritas fam ily (and the corresponding errorbar), as obtained by 
iNesvomv et"ai] j20O3l) . 



Fig.UJa), the random-walkers of the MCMC simulation (triangles) 
practically cover the same region in (Ji,J2) as the real group-A 
members (circles). Moreover, the time evolution of the ratio of the 
standard deviations (j{J2)/ct{Ji), which characterizes the shape of 
the distribution, is reproduced quite well, as shown in Fig.[8lb). 

An additional MCMC simulation with constant (i.e. average 
with respect to and sin Ip) coefficients of diffusion was also per- 
formed. As shown in Fig.[8lb), the value of o{,J2) / (j{Ji) in this 
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Figure 8. Comparison of the time evolution of group A Veritas members 
in (i) a numerical integration of their orbits for 10 Myr, and (ii) an MCMC 
realization, corresponding to the same time interval. Top: the final distribu- 
tion, as taken from the integration (circles) and the MCMC run (triangles). 
Bottom: time evolution of the ratio cr(J2)/o'(Ji). The thick solid (resp. 
dashed) curve corresponds to the MCMC run with variable (resp. constant) 
diffusion coefficients, while the thin, "noisy" one corresponds to the numer- 
ical integration of the equations of motion. 



simulation appears to slowly deviate from the one measured in the 
previous MCMC simulation, as time progresses. However, this de- 
viation is not very large, also compared to the result of the numeri- 
cal integration. Thus, we conclude that an MCMC model with con- 
stant coefficients is adequate for deriving a reasonably accurate es- 
timate of the age of a family, provided that the variations of the local 
diffusion coefficients are not very large. Given this result, we de- 
cided to use average coefficients for the Lixiaohua case. Note also 
that, in the Veritas case, the observed deviation in a{J2) / cr{Ji) 
between the two MCMC models is reflected in the error of r (i.e. 
1.7 Myr vs. 1.2 Myr). 



3.2 The age of the Lixiaohua family 

The Lixiaohua family is another typical outer-belt family, crossed 
by several MMRs. This results into a significant component of fam- 
ily members that follow chaotic trajectories. At the same time, a 
clear 'V -shaped distribution is observed in the (op, H) plane (see 
Fig. |9j, suggesting that the family is old -enough for Yarkovsky to 
have significantly altered its size in ttp. iNesvomv et all ( 120051) in 



this way estimated the age of this family to r = 300 ± 200 
Thus, we choose to study the Lixiaohua family because, on one 
hand, it is relatively old, so Yarkovsky/YORP effects are important, 
but, on the other hand, it has the feature we need (i.e. a significant 
chaotic zone) to test the behaviour of our model on longer time 
scales. Here, we use our MCMC method to derive a more accurate 
estimate of its age, taking into account also the Yarkovsky/YORP 
effects. 

The distribution of the family members in (ap,ep) and 
(ap,sin/p) is shown in Fig. |9] Using velocity cut-off Vc — 
50 m s^^ we find 263 bodies (database as of February 2009) linked 
to the family. The shape of this family is, as in the Veritas case, 
intriguing. For a > 3.15 AU, the family appears to better fit inside 
the equivelocity ellipse shown in the figure, with only a few bodies 
showing a significant excursion in Cp and sin Ip. On the other hand, 
for a < 3.145 AU, the family members occupy a wider area in Bp 
and sin /p. Throughout the family region we find thin, "vertical", 
strips of chaotic bodies, with Lyapunov timefl Tl < 2 x 10* y. 
These strips are associated to different mean motion resonances. 
The most important chaotic domain (hereafter Main Chaotic Zone, 
MCZ) is the one centered around ap ~ 3.146 AU; indicated by 
the grey-shaded area in both plots of Fig.|9] A number of two- and 
three-body MMRs can be associated to the formation of the MCZ, 
such as the 17:8 MMR with Jupiter and the (7, 9, -5) three-body 
MMR. 

Note that the two largest members of this family, (3330) 
Gantrisch and (5900) Jensen, are located just outside the MCZ, as 
indicated by their larger values of Tl (> 2 x lO'* y). In fact, a 
significant group of bodies just outside the MCZ (see Fig. |9ji) has 
higher values of Tl but similar spread in Cp as the MCZ bodies. 
This suggests that bodies around the chaotic zone could have once 
resided therein, evolving towards high/low values of Cp by chaotic 
diffusion. Numerical integrations of the orbits of selected Lixiao- 
hua members for 100 Myr indeed confirmed that bodies could enter 
(or leave) the MCZ. We believe that the distribution of family mem- 
bers on either side of the MCZ is strongly indicative of an interplay 
between Yarkovsky drift in semi-major axis and chaotic diffusion 
in Bp and sin/p, induced by the overlapping resonances; bodies 
can be forced to cross the MCZ, thus receiving a "kick" in Bp and 
sin Ip, before exiting on the other side of the zone. 

A population of ~ 5, 000 fictitious bodies was selected and 
used for calculation of the local diffusion coefficients. Here, we re- 
strict ourselves in calculating coefficients as functions of Op only 
(i.e. averaged in ep and sin Ip). As shown in Fig.[TOl there are sev- 
eral diffusive zones, corresponding to the Iow-Tl strips of Fig. |9] 
However, as in the Veritas case, only one zone appears diffusive 
in both actions; D{J2) is practically zero everywhere outside the 
MCZ (a ~ 3.146 AU), and significant dispersion in proper ele- 
ments is observed only in this zone. 

Given the above results, we conclude that the MCZ family 
members can be used to estimate the age of the family, much like 
the Veritas (5,-2,-2) resonant bodies. Given the fact that random- 
walkers can drift in Op, the way of computing the age is accordingly 
modified. A large number of random walkers, uniformly distributed 
across the whole family region (i.e. the equivelocity ellipses) is 
used. The simulation again stops when 0.3% of MCZ-bodies are 
found to be outside the observed (Ji, J2) borders of the family. 



^ This age was estimated neglecting the initial size of the fa mily. 

^ For details on the computation of Lyapunov times see e.g. lXsiganis et alj 
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Figure 10. The same as Fig. 5, but for the Lixiaohua family. Note that while 
several diffusive peaks are visible in D( J\), only the region of the MCZ 
(Op ~ 3.146 AU) shows significant diffusion in J2. 



3.13 3.135 3.14 3.145 3.15 3.155 3.16 3.165 3.17 
ap [AU] 

Figure 9. Top and middle panels: The same as Fig. 1, but for the Lixi- 
aohua family. The different symbols con'espond to different values of the 
Lyapunov time, Ti^. Triangles correspond to Tj, > 10^ y, empty circles to 
2 X < Tl < 10^ y and solid circles to Tj, < 2 X 10'' y. The equiv- 
elocity curves were obtained for v = AO m'i~^ , f = 80° and lo = 300°. 
Bottom: The distribution of Lixiaohua members in the (ap, H) plane. The 
grey-shaded region denotes the extent of the main chaotic zone (MCZ) in 
ap . The two largest members of the family are also denoted. 



However, the number of MCZ bodies is not constant during the 
simulation, because bodies initially outside (resp. inside) the MCZ 
can enter (resp. leave) that region. Thus, the aforementioned per- 
centage is calculated with respect to the corresponding number of 
MCZ bodies at each time-step. 

The size of the MCZ in the space of proper actions is given by 
cr( Ji) = (8.49±0.51) X 10~* and cr( J2) = (3.17±0.28) x 10"*. 



In order to compute the age of the family we performed 2400 
MCMC runs. This was repeated three times, for three different val- 
ues of thermal conductivity (see above) and once more, neglecting 
the Yarkovsky effect. Given the uncertainty in determining the ini- 
tial size of the family, we repeated the computations for two more 
sets of (5Ji(0), <5J2(0). The results of these computations are pre- 
sented in Fig.[TT] We find an upper limit of 230 Myr for the age 
of the family and a lower limit of ~ 100 Myr. Taking into account 
only the runs performed for our "nominal" initial size and includ- 
ing the Yarkovsky effect, we find the age of Lixiaohua family to be 
155 ± 36 Myr. This value lies towards the lower e nd but comfort- 
ably w ithin the range (300 ± 200 Myr) given by iNesvomv et al.l 
l l2005h . 

We note that, when the Yarkovsky effect is taken into account, 
the age of the family turns out to be longer by ~ 30 Myr (see also 
Fig. lilt . This is a purely dynamical effect, related to the fact that 
bodies can drift towards the MCZ from the adjacent non-diffusive 
regions. However, as more bodies enter the MCZ near its center 
in gp and sin Jp, it takes longer for 0.3% of random walkers to 
diffuse outside the confidence ellipse in (ep,sin/p). At the same 
time, bodies that are initially inside the MCZ can also drift outside, 
to lower/higher values of Op, thus slowing down or even stop diffus- 
ing in Bp and /p. This also explains the large spread in (ep, sin Jp) 
observed for family members located just outside the MCZ. 
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Figure 11. The age of Lixiaohua family, as measured by our MCMC runs. 
The mean value and standard deviation (error-bar) of the age are given, as 
functions of the assumed size of the initial equivelocity ellipse. Four sets of 
points are shown, corresponding to our estimate (i) including Yarkovsky and 
assuming different values of the thermal conductivity (K = 0.01 = dotted 
line, K = 0.1 = dashed line and K = 0.5 = solid thin line), and (ii) 
neglecting Yarkovsky (thick sohd line). 



4 CONCLUSIONS 

We have presented here a refined statistical model for asteroid 
transport, which accounts for the local structure of the phase-space, 
by using variable diffusion coefficients. Also, the model takes into 
account the long-term drift in semi-major axis of asteroids, induced 
by the Yarkovsky/YORP effects. This model can be applied to sim- 
ulate the evolution of asteroid families, also giving rise to an ad- 
vanced version of the "chaotic chronology" method for the deter- 
mination of the age of asteroid families. 

We applied our model to the Veritas family, whose age is well 
constrained from previous works. This allowed us to assess the 
quality and to calibrate our model. We first analyzed the local dif- 
fusion characteristics in the region of Veritas. Our results showed 
that local diffusion coefficients vary by about a factor of 2 — 3 
across the (ep, sin Ip) region covered by the (5,-2,-2) MMR. Thus, 
although local coefficients are needed to accurately model (by the 
MCMC method) the evolution of the distribution of group-A mem- 
bers, average coefficients are enough for a reasonably accurate es- 
timation of the family's age. We note though that the variable co- 
efficients model reduces the error in r by ~ 30%, but requires 
a computationally expensive procedure. Using the variable coeffi- 
cients MCMC model, we found the age of the Veritas family to be 
T = (8.7±1.2) Myr; a res ult in very good agree ment with that of 
iNesvornv et"al] ( l2003h and lTsiganis et al.l ( l2()07h . 

We used our model to estimate also the age of the Lixiao- 
hua family. This family is similar to the Veritas family in many re- 
spects; it is a typical outer-belt family of C-type asteroids, crossed 
by several MMRs. Like the Veritas family, only the main chaotic 
zone (MCZ) shows appreciable diffusion in both eccentricity and 
inclination. On the other hand, this is a much older family and the 
Yarkovsky effect can no longer be ignored. This is evident from the 
distribution of family members, adjacent to the MCZ. Our model 
suggests that the age of this family is between 100 and 230 Myr, 
the best estimate being 155 ± 36 My r. Note that the relativ e error 
is « 23%, i.e. close to the ^ 20% that lTsiganis et al.l j2007h found 
for the Veritas case, using a constant coefficients MCMC model. 

Our model shares some similarities with the Yarkovsky/YORP 



chronology. Both methods are basically statistical and make use 
of the quasi-linear time evolution of certain statistical quantities 
(either the spread in Aa or the dispersion in and sin/p), de- 
scribing a family. There are, however, important differences. The 
Yarkovsky/YORP chronology method works better for older fami- 
lies and the age estimates are more accurate for this class of asteroid 
families (provided there are no other important effects on that time 
scale). On the other hand, for our method to be efficient, we need 
that diffusion is fast enough to cause measurable effects, but slow 
enough so that most of the family members are still forming a ro- 
bust family structure (i.e. there is no dynamical "sink" that would 
lead to a severe depletion of the chaotic zone). Thus, our model 
can be applied to a limited number of families that reside in com- 
plex phase-space regions, but, in the same time, this is the only 
model that takes into account the chaotic dispersion of these fam- 
ilies. There are at least a few families for which both chronology 
methods can be applied, thus leading to more reliable age estimates, 
as well as to a direct comparison of the two different chronologies. 
For example, the families of (20) Massalia and (778) Theobalda 
would be good test cases. We, however, reserve this for future work. 

An important advantage of the model is that it can be used to 
estimate the physical properties of a dynamically complex aster- 
oid family, provided that its ag e is known by indepen dent means 
(e.g. by applying the method of iNesvomv et al.l j2003h to the reg- 
ular members of the family). A large number of MCMC runs can 
be performed at low computational cost, thus allowing a thorough 
analysis of the physical parameters of family members or the prop- 
erties of the original ejection velocities field that better reproduce 
the currently observed shape of the family. 
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